#!/usr/bin/env python
#########################################################################################
#
# PTCLMmkdata
#
# Python script to create datasets to run point simulations of CLM5
# using Tower Datasets for Ameriflux tower sites, using the CESM2
# framework.
#
# Python script originally created by: 
#
#  Daniel M. Riccciuto, Dali Wang, Peter E. Thornton, Wilfred M. Poist
#  of Environmental Sciences Division, Oak Ridge National Lab.
#
#  Quinn Thomas
#  of Cornell University
#
#  Modified by Erik Kluzek (NCAR) to be incorporated as a standard part of CLM5.
#
#  For help on PTCLMmkdata type:
#
#   PTCLMmkdata --help
#
#  Also see the README file
#
#  Requirements:
#
#  python, UNIX shell, NCL (NCAR Command Language),
#  GNU make, Fortran compiler, C compiler
#
# NOTE:  mksurfdata_map, and gen_domain  must be compiled!
#           You should only have to compile them once.
#           you must also have ncl installed.
#
#########################################################################################
description = 'Python script to create datasets to run single point simulations with tower site data.'
import os, csv, time, re, sys, shlex, subprocess
from   xml.sax.handler import ContentHandler
from   xml.sax         import make_parser 

######  THE ERROR FUNCTION
##############################################################

def error( desc ):
     "error function"
     print( "ERROR("+sys.argv[0]+"):: "+desc )
     os.abort()

######  SET SOME VARIABLES ##############################################################

#configure case options          
#run time defaults
defSitesGroup = "PTCLMDATA" #default site group name

global ptclm_dir
stdout    = os.popen("pwd")
cwd       = os.path.abspath( stdout.read().rstrip( ) )
dirname   = os.path.dirname(sys.argv[0])
if ( dirname == "" ):
   ptclm_dir = cwd
else:
   ptclm_dir = os.path.abspath(dirname)

base_ctsm     = os.path.abspath(ptclm_dir+"/../.." )
wrkdir        = cwd
cesm_input    = " "
filen         = " "
mydatadir     = ptclm_dir+"/mydatafiles"
clmnmlusecase = "2000_control"
stdout        = os.popen( "date +%y%m%d" );
sdate         = stdout.read().rstrip( );

######  GET VERSION INFORMATION #########################################################

if sys.version_info < (2, 4):
    error( "The version of Python being used is too old for PTCLMmkdata" )


tagvers = ""
clog = open( ptclm_dir+"/ChangeLog", "r" );
for line in clog:
   if ( line.find("Tag: ",0) == 0 ): 
      n = line.count("")
      tagvers = line[5:n-2]
      break

clog.close
version=tagvers

### PARSE THE COMMAND LINE INPUT ########################################################

from optparse import OptionParser, OptionGroup

#parse arguments
cmdline = ""
for arg in sys.argv:
    cmdline = cmdline+arg+" "
parser = OptionParser( usage="%prog [options] -d inputdatadir -s sitename", description=description, version=version )
required = OptionGroup( parser, "Required Options" )
required.add_option("-d", "--cesmdata", dest="cesm_input", default=" ", \
                  help="Location of CCSM input data")
required.add_option("-s", "--site", dest="mysite", default="none", \
                  help="Site-code to run, FLUXNET code (-s list to list valid names)")
parser.add_option_group(required)
options  = OptionGroup( parser, "Configure and Run Options" )
options.add_option("--crop", dest="crop", \
                  action="store_true", help = \
                  "Create datasets and run with prognostic crop on")
options.add_option("--no-crop", dest="crop", \
                  action="store_false", default=False, help = \
                  "Create datasets and run without prognostic crop on (DEFAULT)")
options.add_option("--ctsm_root", dest="base_ctsm", \
                  default=base_ctsm, help = \
                  "Root CTSM directory (top level directory CTSM src bld doc cime_config subdirs)")
options.add_option("--cime_root", dest="base_cime", \
                  default=" ", help = \
                  "Root CIME directory (top level directory of CIME infrastructure with " \
                  + " scripts, doc, utils, tools, src, config subdirs)")
options.add_option("--debug", dest="debug", action="store_true", default=False, \
                  help="Flag to turn on debug mode so won't run, but display what would happen")
options.add_option("--sdate", dest="sdate", default=sdate, \
                  help="Use entered date string in all files"+\
                       " (use the given date string in place of the current date:"+sdate+")" )
options.add_option("--clmnmlusecase", dest="clmnmlusecase", default=clmnmlusecase, \
                  help="CTSM namelist use case to use (default:"+clmnmlusecase+")" )
options.add_option("--list", dest="list", default=False, action="store_true", \
                  help="List all valid: sites")
options.add_option("--mydatadir", dest="mydatadir", default=mydatadir \
                  ,help="Directory of where to put your data files (files will be under subdirectories for each site)"+\
                        " (default: "+mydatadir+")" )
options.add_option("--donot_use_tower_yrs",action="store_false",\
                  dest="use_tower_yrs",default=True,\
                  help="Do NOT use the data years that correspond to the tower years "+\
                       "(when you plan on using global forcing)" )
options.add_option("--quiet", action="store_true", \
                  dest="quiet", default=False, \
                  help="Print minimul information on what the script is doing")
options.add_option("--cycle_forcing", action="store_true", \
                  dest="cycle_forcing", default=False, \
                  help="Cycle over the forcing data rather than do one run through (modifies start/end year to get this to work)")
options.add_option("--verbose", action="store_true", \
                  dest="verbose", default=False, \
                  help="Print out extra information on what the script is doing")
parser.add_option_group(options)

indatgengroup = OptionGroup( parser, "Input data generation options", \
                  "These are options having to do with generation of input datasets.  " )
parser.add_option_group(indatgengroup)
indatgengroup.add_option("--pftgrid", dest="pftgrid", help = \
                  "Use pft information from global gridded file (rather than site data)", \
                  action="store_true", default=False)
indatgengroup.add_option("--soilgrid", dest="soilgrid", help = \
                  "Use soil information from global gridded file (rather than site data)",\
                   action="store_true", default=False)
indatgengroup.add_option("--map_gdate", dest="map_gdate", default=sdate, \
                  help="Use existing mapping files with the given date string rather than create new ones with current date"+\
                       " (if mapping files do NOT exist with this date, the script will abort)" )
indatgengroup.add_option("--mksurfdata_opts", dest="mksurfdata_opts", help = \
                  "Options to send directly to mksurfdata_map",\
                   default="")
versiongroup  = OptionGroup( parser, tagvers )
parser.add_option_group(versiongroup)

(options, args) = parser.parse_args()
if len(args) != 0:
    parser.error("incorrect number of arguments")

### END PARSE THE COMMAND LINE INPUT ####################################################

### SOME FUNCTIONS    ###################################################################

def system( cmd ):
     "system function with error checking and debug prining"
     global plev

     if plev>0: print( "Run command: "+cmd )
 
     # Check if this is a command to always do regardless of debug
     cmdsallow = [ "create_newcase", "mkdir", "mv", "cat", "which", "chmod", "touch", "mksurfdata.pl" ]
     allowed_cmd = False
     for allow_cmd in cmdsallow:
        if ( cmd.find( allow_cmd ) > 0 ):
          allowed_cmd = True

     # Error check that command exists
     if ( not options.debug or allowed_cmd ):
        firstspace = cmd.index(" ");
        if ( firstspace == -1 ):
           justcmd = cmd
        else:
           justcmd = cmd[:firstspace]

        if ( cmd.index("/") != -1 ):
           if ( not os.path.exists(justcmd) ): 
              error( "Error command does NOT exist: "+justcmd );
        else:
           rcode = os.system( "which "+justcmd )
           if ( rcode != 0 ):
              error( "Error command is NOT in path: "+justcmd )

     # Now actually run the command
     if ( not options.debug or allowed_cmd ):
        if ( options.debug and cmd.find( "mksurfdata.pl" ) > 0):
           rcode = os.system( cmd+" --debug --allownofile" )
        else:
           rcode = os.system( cmd )
     else: rcode = 0
     if ( rcode != 0 ):
        error( "Error running command: "+cmd )
        if ( os.path.isfile(filen) ):
           output = open( filen,'a')
           output.write(cmd+"\n")
           output.close

def queryFilename( queryopts, filetype ):
    "query the XML database to get a filename"
    global abs_base_ctsm
    query = abs_base_ctsm+"/bld/queryDefaultNamelist.pl -silent " \
             +"-justvalue "
    if ( cesm_input != " " ): 
       query = query + " -csmdata "+cesm_input
    cmd = query+queryopts+" -var "+filetype
    file = os.popen( cmd )
    filename = file.read() 
    if ( file.close() != None ):
       print( "Query = "+cmd )
       error( "Error getting file from XML database" )
    # Remove the trailing new line from the filename
    if ( (filename == None) or (filename == "") ): 
       print( "Query = "+cmd )
       error( "Trouble finding file from XML database: "+filetype )
    return( filename.replace( "\n", "" ) )

def setup_case_files( ):
     "Setup the user_nl_clm and shell_commands files"

     global data_dir
     global cmdline
     filex = data_dir+"/shell_commands"
     output = open( filex,'w')
     output.write("# shell commands to execute xmlchange commands written by PTCLMmkdata:\n")
     output.write("# "+cmdline+"\n")
     output.close
     system( "/bin/chmod +x "+filex )
     usernlclm = data_dir+"/user_nl_clm"
     output = open( usernlclm,'w')
     output.write("! user_nl_clm namelist options written by PTCLMmkdata:\n")
     output.write("! "+cmdline+"\n")
     output.close
     return( filex, usernlclm )

def xmlchange_env_value( filex, var, value, note="", append=False ):
     'Function to set the value of a variable in one of the env_*.xml files'
     change = "./xmlchange"
     if ( append ):
        change = change + " --append"
     cmd = change+" "+var+"="+value
     output = open( filex,'a')
     if ( note != "" ):
        output.write("# "+note+"\n")
     output.write(cmd+"\n")
     output.close

def write_datm_namelistdefaults_file( dir ):
     "Write namelist_defaults_datm.xml file"
     global data_dir
     datm_src_dir = data_dir+"/SourceMods/src.datm";
     os.system( "/bin/mkdir -p "+datm_src_dir )
     file = datm_src_dir+"/namelist_defaults_datm.xml"
     output = open( file,'w')
     filestrings = ( '<?xml version="1.0"?>', \
     ' ', \
     '<?xml-stylesheet type="text/xsl" href="namelist_defaults.xsl"?>', \
     ' ', \
     '<namelist_defaults>', \
     ' ', \
     '</namelist_defaults>' )
     # write out file
     for line in filestrings:
        output.write(line+"\n")
        # Add data directory for files
        if ( line.find( "<namelist_defaults>" ) != -1 ):
           output.write('\n')
           value = '<strm_datdir stream="CLM1PT.CLM_USRDAT">'+dir+'</strm_datdir>\n'
           output.write(value)

     output.close

def find_filename_created( wildcard, desc ):
     "Find the filename of the file that was just created"

     global options
     if ( not options.debug ):
        # If NOT debug mode, get the filename from a directory listing
        stdout   = os.popen( "ls -1t1 "+wildcard+" | head -1" );
        filename = stdout.read().rstrip( );
        if ( not os.path.exists( filename) ): error( "filename does NOT exist:"+wildcard )
     else:
        # For debug mode, create a file with current date replacing any wildcards
        filename = wildcard.replace( "*", options.sdate )
        os.system( "touch "+filename )

     print( desc+" = "+filename )
     return( filename )


if sys.version_info < (2, 5):
   def rpartition( string, sep ):
       'Reverse order of dividing string by seperator'
       before = string[0:string.rfind(sep)];
       after  = string[before.count(""):];
       return ( before, sep, after )

###### SET OPTIONS BASED ON INPUT FROM PARSER  ##########################################

mysite     = options.mysite
SitesGroup = defSitesGroup
infohelp   = "\n\n Use --help option for help on usage.\n";
if(options.list):
    mysite    = "list"
if ( mysite == "none" ): parser.error("sitename is a required argument, set it to a valid value"+infohelp )
if ( options.verbose and options.quiet ):
    parser.error( "options quiet and verbose are mutually exclusive"+infohelp )

if (   options.verbose ): plev = 2
elif ( options.quiet   ): plev = 0
else:                     plev = 1

sitedata=SitesGroup+"_sitedata.txt"
soildata=SitesGroup+"_soildata.txt"
pftdata=SitesGroup+"_pftdata.txt"


if plev>0: print( "---------------- PTCLMmkdata version "+str(version)+"-----------------------------\n" )
if plev>0: print( "   "+cmdline+"\n" )
if plev>0: print( "   OPTIONS:\n" )
if plev>0: print( "Site name:\t\t\t\t\t\t"+mysite+"\n" )

base_ctsm = options.base_ctsm
if base_ctsm == " ":
    #assume base directory is three levels up from where PTCLM script
    #  is executed, if not specified
    stdout    = os.popen("cd "+ptclm_dir+"/../../..; pwd")
    base_ctsm = os.path.abspath( stdout.read().rstrip( ) )

abs_base_ctsm = os.path.abspath( base_ctsm )
if plev>0: print( "Root CTSM directory:\t\t\t\t\t"+abs_base_ctsm )

base_cime = options.base_cime
if base_cime == " ":
    #assume base directory is tunder the CTSM directory
    #  is executed, if not specified
    stdout    = os.popen("cd "+abs_base_ctsm+"/cime; pwd")
    base_cime = os.path.abspath( stdout.read().rstrip( ) )

abs_base_cime = os.path.abspath( base_cime )
if plev>0: print( "Root CIME directory:\t\t\t\t\t"+abs_base_cime )

if plev>0: print( "** Surface data file will be built using site-level data " + 
          "when available unless otherwise specified ** \n" )
if plev>0: print( "\tExtract PFT data from gridded files:\t\t"+str(options.pftgrid) )
if plev>0: print( "\tExtract soil data from gridded files:\t\t"+str(options.soilgrid) )

###### END SET OPTIONS BASED ON INPUT FROM PARSER  ######################################

########## GET SITE LAT, LON, AND TOWER MET YEARS #######################################

siteDir = ptclm_dir+"/"+"PTCLM_sitedata"
#get lat/lon, start/end years from sitedata file
if plev>0: print( "\nOpen Site data file: "+siteDir+"/"+sitedata+"\n" )
sitepath = siteDir+"/"+sitedata
AFdatareader = csv.reader(open(sitepath, "rb"))
if ( mysite == "list" ):  plev = 2
found=False
for row in AFdatareader:
    if plev>1: print( " site = %9s name: %-55s Region: %12s Campaign: %s" % ( row[0], row[1], row[2], row[10] ) )
    if row[0] == mysite:
        found=True
        lon=float(row[3])
        if (lon < 0):
            lon=360.0+float(row[3]) 
        lat=float(row[4])
        startyear=int(row[6])
        endyear=int(row[7])
        alignyear = int(row[8])
        timestep  = int(row[9])

# Exit early for list options
if ( mysite == "list" ): 
  exit()
if ( not found ):
  parser.error( "Entered site is NOT in the list of valid sites: "+mysite )

# inputdata directory -- set after list options
cesm_input=options.cesm_input
if cesm_input == " ":
   parser.error( "inputdatadir is a required argument, set it to the directory where you have your inputdata"+infohelp )
if plev>0: print( "CESM input data directory:\t\t\t\t"+cesm_input )
#define data and utility directories
mask          = "navy"
clmusrdatname = "1x1pt_"+mysite
clmusrdat     = " -usrname "+clmusrdatname
clmres        = clmusrdatname
clmmask       = "navy"
myres         = "CLM_USRDAT"  #single-point mode (don't change)

clm_tools   = abs_base_ctsm+'/tools'
if ( not os.path.exists( clm_tools ) ): 
   error( "clm tools directory does NOT exist: "+clm_tools )
gen_dom_dir = abs_base_cime+'/tools/mapping/gen_domain_files'
if ( not os.path.exists( gen_dom_dir ) ): 
   error( "generate domain directory does NOT exist: "+gen_dom_dir )
mkmapgrd_dir= clm_tools+'/mkmapgrids'
if ( not os.path.exists( mkmapgrd_dir ) ): 
   error( "make map grid directory does NOT exist: "+mkmapgrd_dir )
mkmapdat_dir= clm_tools+'/mkmapdata'
if ( not os.path.exists( mkmapdat_dir ) ): 
   error( "make map data directory does NOT exist: "+mkmapdat_dir )
clm_input   = cesm_input+'/lnd/clm2'
datm_input  = cesm_input+'/atm/datm7'

mydata_dir  = os.path.abspath( options.mydatadir )
data_dir    = mydata_dir+"/"+clmusrdatname
if ( not os.path.exists( data_dir ) ): os.system( "/bin/mkdir -p "+data_dir )

if plev>0: print( "----------------------------------------------------------------\n" )

############# WRITE OUT README FILE ON DATA ##############################################

opt = " "

clmnmlusecase    = options.clmnmlusecase

filen = data_dir+"/README.PTCLM"
if plev>0: print( "Write "+filen+" with command line" )
output = open( filen,'w')
output.write(cmdline+"\n")
output.close

############# GET SIM_YEAR, RCP and SIM_YEAR_RANGE based on USE-CASE ####################
############# CLM configure ensures naming conventions are followed  ####################
############# And setup Query options based on them #####################################

if (   clmnmlusecase.endswith("_transient") ):
     transient = re.search('^([0-9]+-[0-9]+)_*(.*)_(transient$)',   clmnmlusecase )
     if ( transient ):
        sim_year_range = transient.group(1)
        sim_year       = re.search( '^([0-9]+)-',    transient.group(1) ).group(1)
     elif ( clmnmlusecase.startswith("20thC_") ):
        sim_year_range = "1850-2000"
        sim_year       = "1850"
     else:
        error( "Can not parse use-case name, does not follow conventions: "+clmnmlusecase )

     if ( sim_year_range == "1850-2000" ): actual_sim_year_range = "1850-2015"
     else:                                 actual_sim_year_range = sim_year_range
elif ( clmnmlusecase.endswith("_control") ):
     control        = re.search( '^([0-9]+)_', clmnmlusecase )
     if ( not control ):      error( "Can NOT parse use-case name does NOT follow conventions: "+clmnmlusecase )
     sim_year       = control.group(1)
     if ( sim_year == None ): error( "Trouble finding sim_year from:"+clmnmlusecase )
     sim_year       = str(sim_year)
     sim_year_range = "constant"
elif ( clmnmlusecase.endswith("_pd") or clmnmlusecase == "UNSET" ):
     sim_year       = "2000"   
     sim_year_range = "constant"
else:
     error( "Can not parse use-case name:, does not follow conventions: "+clmnmlusecase )

landuse_timeseries_type = "hist"

qoptionsbase   = " -options mask="+mask
   
qoptions       = qoptionsbase+",sim_year="+sim_year+",sim_year_range="+sim_year_range;
queryOpts      = " -onlyfiles -res "+clmres+clmusrdat+qoptions
queryOptsNousr = qoptions
queryOptsNavy  = " -res 0.33x0.33 "+qoptions

#
# If you are trying to cycle the forcing years you need to be careful about
# the number of years cycling over and taking leap years into account.
#
if ( options.cycle_forcing ):
    numyears = endyear - startyear + 1
    numfour = int(numyears/4)
    # If have three years or less (numfour = 0) just repeat first year  
    # unless first year is leap year then use next year.
    # Since just using one year that is not a leap year endyear is startyear
    if (numfour == 0):
      if (startyear % 4 == 0):
        startyear = startyear + 1

      endyear  = startyear
    else:
      endyear = startyear + numfour * 4 - 1

    # Use alignyear from file for cycle_forcing case
else:
    # When NOT cycling forcing, use start year for the align year
    alignyear = startyear

if (options.crop == True):
   mkcrop = " -crop"
else:
   mkcrop = " -no-crop"

####### ANY OTHER LAST SETTINGS BEFORE CREATING DATASETS ################################

#####  ENV XML CHANGES ##################################################################
filex, usernlclm = setup_case_files( )

if ( clmusrdatname != "" ):
   xmlchange_env_value( filex, "CLM_USRDAT_NAME", clmusrdatname )

if(options.use_tower_yrs):
    xmlchange_env_value( filex, "DATM_CLMNCEP_YR_START", str(startyear) )
    xmlchange_env_value( filex, "DATM_CLMNCEP_YR_END",   str(endyear) )

xmlchange_env_value( filex, "MPILIB", "mpi-serial", note="Comment this out if NINST_LND is greater than 1 (see: http://bugs.cgd.ucar.edu/show_bug.cgi?id=2521)" )

############# BEGIN CREATE POINT DATASETS ###############################################


if plev>0: print("Making input files for the point (this may take a while if creating transient datasets)")

os.chdir(data_dir)
#make map grid file and atm to ocean map ############################################
if plev>0: print( "Creating map file for a point with no ocean" )
print( "lat="+str(lat) )
ptstr = str(lat)+","+str(lon)
if ( os.system( "which ncl" ) != 0 ): error( "ncl is NOT in path" )  # check for ncl
system(mkmapdat_dir+"/mknoocnmap.pl -p "+ptstr+" -name "+clmres+" > "+data_dir+"/mknoocnmap.log")
mapfile = find_filename_created( data_dir+"/map_"+clmres+"_noocean_to_"+clmres+"_"+"nomask_aave_da_*.nc", "mapfile" )
scripgridfile  = find_filename_created( data_dir+"/SCRIPgrid_"+clmres+"_nomask_c*.nc", "scripgridfile" )

#make domain file needed by datm ####################################################
if plev>0: print( "Creating data domain" )
cmd = gen_dom_dir+"/gen_domain -m "+mapfile+" -o "+clmmask+" -l "+clmres+" -c 'Running gen_domain from PTCLMmkdata' > "+data_dir+"/gen_domain.log"
system(cmd);
domainfile  = find_filename_created( "domain.lnd."+clmres+"_"+clmmask+".*.nc", "domainfile" )

#make surface data and dynpft #######################################################
if plev>0: print( "\n\nRe-create surface dataset:\t" )
if ( sim_year_range == "constant" ):
   mksrfyears = sim_year
else:
   mksrfyears = sim_year_range

#make mapping files needed for mksurfdata_map #######################################

mapdir = data_dir
if ( options.map_gdate == options.sdate ):
  # mkmapdata.sh remembers where it is (although it starts over for a new date)
  if plev>0: print( "\n\nRe-create mapping files for surface dataset:" )
  cmd = "export REGRID_PROC=1; "+mkmapdat_dir+"/mkmapdata.sh --gridfile "+scripgridfile+" --res "+clmres+" --gridtype regional -v > "+mapdir+"/mkmapdata.log";
  system(cmd);
else:
  mksrfmapfile  = find_filename_created( mapdir+"/map_*"+"_c"+options.map_gdate+".nc", "mksrfmapfile" )
  if ( not os.path.exists( mksrfmapfile ) ): error( "mapping files with gdate of "+ \
       options.map_gdate+" do NOT exist, bad value for --map_gdate option"     )

# --- use site-level data for mksurfdata_map when available ----
#PFT information for the site
if (options.pftgrid == False):
    if plev>0: print( "Replacing PFT information in surface data file" )
    os.chdir(siteDir)
    AFdatareader = csv.reader(open(pftdata, "rb"))
    os.chdir(data_dir)
    pft_frac=[0,0,0,0,0]
    pft_code=[0,0,0,0,0]
    found=0
    for row in AFdatareader:
        if plev>1: print( " site = %9s" % row[0] )
        if row[0] == mysite:
            found=1
            output=open("./tempsitePFT.txt","w")      
            output.write(' '.join(row[1:11]))
            output.close()
            for thispft in range(0,5):
                pft_frac[thispft]=float(row[1+2*thispft])
                pft_code[thispft]=int(row[2+2*thispft])
    if ( found == 0 ):
       error( "Did NOT find input sitename:"+mysite+" in pftdata:"+pftdata+ \
              " run with pftgrid instead")
    # Find index of first zero
    for i in range(0,len(pft_frac)):
       if ( pft_frac[i] == 0.0 ):
          nzero = i
          break
    pftopts=" -pft_frc \""+str(pft_frac[0:nzero])+'"' \
               " -pft_idx \""+str(pft_code[0:nzero]) +'"' + mkcrop
else: 
    pftopts=""
     
#Read in the soil conditions for the site #######################################
if (options.soilgrid == False):

    #soil information
    os.chdir(siteDir)
    if plev>0: print( "Replacing soil information in surface data file" )
    AFdatareader = csv.reader(open(soildata, "rb"))
    os.chdir(data_dir)
    found=0
    for row in AFdatareader:
        if plev>1: print( " site = %9s" % row[0] )
        if row[0] == mysite:
            found=1
            output=open("./tempsitesoil.txt","w")
            output.write(' '.join(row[1:7]))
            output.close()
            # The first three items are NOT used
            soil_depth = float(row[1])  # This is ignored
            n_layers   = int(row[2])    # This is ignored
            layer_depth = float(row[3]) # This is ignored
            sandpct     = float(row[4])
            claypct     = float(row[5])
    if ( found == 0 ):
       error( "Did NOT find input sitename:"+mysite+" in soildata:"+soildata+ \
              " run with soilgrid instead")
    if plev>0: print( " sandpct="+str(sandpct)+" claypct="+str(claypct) )
    soilopts=" -soil_cly "+str(claypct)+" -soil_snd "+str(sandpct)
else: soilopts=""
#----- create dynamic pft input file --------------- ############################
if (options.pftgrid == False) and (sim_year_range != "constant"):

    if plev>0: print( "Creating site-specific dynamics PFTs and harvesting" )

    landuse_timeseries_site_filename = siteDir + \
                           mysite + "_dynpftdata.txt"

    # only set dynpft file if the file exists
    if ( os.path.exists( landuse_timeseries_site_filename ) ):
       if plev>0: print( "Transition PFT file exists, so using it for changes in PFT" )
           # Convert the file from transition years format to mksurfdata_map landuse_timeseries_ format
       cnv = siteDir + \
             "/cnvrt_trnsyrs2_landuse_timeseries_txtfile.pl " + \
             landuse_timeseries_site_filename+" "+sim_year_range
       landuse_timeseries_outfile = data_dir+"/landuse_timeseries_"+mysite+".txt"
       system( cnv+" > "+landuse_timeseries_outfile )
       dynpftopts = " -dynpft "+landuse_timeseries_outfile
    else:
       error( "Transition PFT file does NOT exist for this site, create one, use --pftgrid, or choose a non transient use-case" )
           
else: 
    dynpftopts = ""

# Now run mksurfdata_map  ###########################################################
mksurfopts = "-res usrspec -usr_gname "+clmres+" -usr_gdate "+options.map_gdate+ \
             " -usr_mapdir "+mapdir+" -dinlc "+cesm_input+" -y "+mksrfyears+ \
             soilopts+pftopts+dynpftopts+" "+options.mksurfdata_opts
system(clm_tools+"/mksurfdata_map/mksurfdata.pl "+mksurfopts+" > "+data_dir+"/mksurfdata_map.log")

surffile  = find_filename_created( data_dir+"/surfdata_"+clmres+"*_simyr"+sim_year+"_*.nc",  "surface file"     )
logfile   = find_filename_created( data_dir+"/surfdata_"+clmres+"*_simyr"+sim_year+"_*.log", "surface log file" )
if ( sim_year_range != "constant" ):
   landuse_timeseries_file = find_filename_created( data_dir+"/landuse.timeseries_"+clmres+"_"+landuse_timeseries_type+"*_simyr"+actual_sim_year_range+"_*.nc", "landuse_timeseries_file" )
# rename files with clm version in the filename
mkopts = ""
if (options.pftgrid         == True):  mkopts += "_pftgrd"
if (options.soilgrid        == True):  mkopts += "_soigrd"
if (options.mksurfdata_opts != ""  ):  mkopts += "_"+options.mksurfdata_opts.replace(" ","+")


   
####### END CREATE POINT DATASETS #######################################################


###### SET ENV_RUN.XML VALUES ###########################################################
   
os.chdir(data_dir)
xmlchange_env_value( filex, "ATM_DOMAIN_PATH", data_dir   )
xmlchange_env_value( filex, "LND_DOMAIN_PATH", data_dir   )
xmlchange_env_value( filex, "ATM_DOMAIN_FILE", domainfile )
xmlchange_env_value( filex, "LND_DOMAIN_FILE", domainfile )
xmlchange_env_value( filex, "CLM_BLDNML_OPTS", "'-mask "+mask+mkcrop+"'", append=True );

xmlchange_env_value( filex, "CALENDAR",        "GREGORIAN" )
xmlchange_env_value( filex, "DOUT_S",          "FALSE" )
hist_nhtfrq = 0
hist_mfilt  = 1200

atm_ncpl = int((60 // timestep) * 24)
xmlchange_env_value( filex,    "ATM_NCPL", str(atm_ncpl) )
if(options.use_tower_yrs):
   xmlchange_env_value( filex, "RUN_STARTDATE", str(alignyear)+"-01-01" )
   xmlchange_env_value( filex, "DATM_CLMNCEP_YR_ALIGN", str(alignyear) )

xmlchange_env_value( filex, "DIN_LOC_ROOT",         cesm_input )
xmlchange_env_value( filex, "DIN_LOC_ROOT_CLMFORC", mydata_dir )

####  NAMELIST DEFAULTS FILE MODIFICATIONS ##############################################
datm_dir = data_dir+"/CLM1PT_data"
if ( os.path.isdir(datm_dir) ):
   write_datm_namelistdefaults_file( datm_dir )

####  SET NAMELIST OPTIONS ##############################################################
output = open(usernlclm,'a')
output.write(   " fsurdat = '"+surffile+"'\n" )
if (sim_year_range != "constant"):
   output.write(   " flanduse_timeseries = "+landuse_timeseries_file+"\n" )
output.write(   " hist_nhtfrq = "+str(hist_nhtfrq)+"\n" )
output.write(   " hist_mfilt  = "+str(hist_mfilt)+"\n" )
output.close()
if plev>1: os.system( "/bin/cat user_nl_clm" )

###### END SET Spinup and ENV_RUN.XML VALUES ############################################

if plev>0: print( "Data created successfully in "+data_dir+"\n" )

###   END PTCLM SCRIPT ####################################################################

